<p>
  Implement a 2D max pooling operation for image/feature map downsampling.
  The program should take an input tensor and produce an output tensor by applying max pooling with specified kernel size, stride, and padding.
</p>

<h2>Implementation Requirements</h2>
<ul>
  <li>External libraries are not permitted</li>
  <li>The <code>solve</code> function signature must remain unchanged</li>
  <li>The final result must be stored in tensor <code>output</code></li>
</ul>

<h2>Max Pooling Operation</h2>
<p>
  For each output position (n, c, h_out, w_out), compute the maximum value over the corresponding input window:
  <br>
  <code>output[n, c, h_out, w_out] = max(input[n, c, h:h+kernel_size, w:w+kernel_size])</code>
  <br>
  where h = h_out * stride and w = w_out * stride
</p>

<h2>Example 1:</h2>
<pre>
Input:  input = [[[[1.0, 2.0, 3.0],
                   [4.0, 5.0, 6.0],
                   [7.0, 8.0, 9.0]]]]
        kernel_size = 2
        stride = 1
        padding = 0
Output: output = [[[[5.0, 6.0],
                    [8.0, 9.0]]]]
</pre>

<h2>Example 2:</h2>
<pre>
Input:  input = [[[[1.0, 2.0, 3.0, 4.0, 5.0],
                   [6.0, 7.0, 8.0, 9.0, 10.0],
                   [11.0, 12.0, 13.0, 14.0, 15.0],
                   [16.0, 17.0, 18.0, 19.0, 20.0],
                   [21.0, 22.0, 23.0, 24.0, 25.0]]]]
        kernel_size = 3
        stride = 1
        padding = 1
Output: output = [[[[13.0, 14.0, 15.0],
                    [18.0, 19.0, 20.0],
                    [23.0, 24.0, 25.0]]]]
</pre>

<h2>Constraints</h2>
<ul>
  <li>1 ≤ N ≤ 100 (batch size)</li>
  <li>1 ≤ C ≤ 512 (channels)</li>
  <li>1 ≤ H, W ≤ 1024 (height, width)</li>
  <li>1 ≤ kernel_size ≤ 16</li>
  <li>1 ≤ stride ≤ 16</li>
  <li>0 ≤ padding ≤ 16</li>
  <li>Input and output tensors use float32 precision</li>
</ul>